#include <stdlib.h>
#include "alg/random.h"
#include "alg/rmat.h"


void rmat_edge (int64_t * iout, int64_t * jout, int SCALE, double A, double B, double C, double D, dxor128_env_t * env)
{
    int64_t i = 0, j = 0;
    int64_t bit = ((int64_t) 1) << (SCALE - 1);

    while (1) {
        const double r = dxor128(env);
        if (r > A) {                /* outside quadrant 1 */
            if (r <= A + B)           /* in quadrant 2 */
                j |= bit;
            else if (r <= A + B + C)  /* in quadrant 3 */
                i |= bit;
            else {                    /* in quadrant 4 */
                j |= bit;
                i |= bit;
            }
        }
        if (1 == bit)
            break;

        /*
        Assuming R is in (0, 1), 0.95 + 0.1 * R is in (0.95, 1.05).
        So the new probabilities are *not* the old +/- 10% but
        instead the old +/- 5%.
        */
        A *= (9.5 + dxor128(env)) / 10;
        B *= (9.5 + dxor128(env)) / 10;
        C *= (9.5 + dxor128(env)) / 10;
        D *= (9.5 + dxor128(env)) / 10;
        /* Used 5 random numbers. */

        {
            const double norm = 1.0 / (A + B + C + D);
            A *= norm;
            B *= norm;
            C *= norm;
        }
        /* So long as +/- are monotonic, ensure a+b+c+d <= 1.0 */
        D = 1.0 - (A + B + C);

        bit >>= 1;
    }

    /* Iterates SCALE times. */
    *iout = i;
    *jout = j;
}

